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See the Figure on the cover page. 


A regular pentagon 
forms a sort of arrow head: 


A chain of pentagons is formed by 

appending a new pentagon to an 

old one on one of the two sides of the 

arrow head. The choice of side is determined : 

by the odd or even nature of the digits in the decimal | 
The first 29 links of the dragon are } 

shown in the cover Figure. : 


Where is the center of the 1000th link of 
the dragon? Assume that the START pentagon is 
centered at the origin and that the pentagons have 
unit sides. 


The drawing on page 2 shows the information needed to 


At the stage where pentagon E is known (centered at : 
C,D and pointing in direction Q), the calculation must : 
determine the coordinates A,B of pentagon F, which : 
points in direction R, 54° clockwise from Q. The 

case shown is for a turn to the right, as dictated 

by an even digit ina. 


The 1000 digits of wshown at the left were keypunched 
50 digits per card, in columns 1-50. The listing 
shows 55 lines of 18 digits and one line of 10 digits. 


(1) Write a program to read the 20 cards and print 
the 56 lines. The program can be in assembly 
language, or any higher level language. 


(2) Indicate, for any such program, how to change 
the parameter 18 to any other value. 


(3) Indicate a test procedure for the program. 
PROBLEM 149 
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THE DISTRIBUTION OF NUMBERS--COMPUTER THEORY 


by R. W. Hamming 


Bell Laboratories Murray Hill, New Jersey 


It is reasonably well known that the fraction of 
i physical constants whose mantissas (in scientific notation) 
i are less than x is 


* Jogex 
D(x) = log b 


i where b is the base of the number system, usually 10. 

i Thus D(x) is the probability of observing a mantissa which 
i is less than x when picking a "random physical constant." 

i A theoretical derivation for this observed cumulative 

: probability distribution can be made based on the 

+ assumption that the distribution is not changed when the 

: individual members are multiplied by some non-zero constant 
ik; that is, 


D(kx) = D(k) + D(x). 


: This same distribution arises when we examine the i @ 
: question of the distribution of floating point numbers 
: in a computing machine. In particular, we shall show 
: that the processes of multiplication and division tend 
i to produce this distribution. 


: As a first step we consider a one decimal digit 
? computer, and examine the 81 possible products rounded to 
i one digit (of the products ending in 5 we round up half 

: the time and down half the time, leaving the single 

: product 5x5 to round up). From this experiment, we get 
i the following table of products: 
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If we regard each of the original numbers as being 
equilikely, then the 81 products are all equilikely, and 
column 2 of the next table gives the corresponding 
frequencies, while column 3 gives the cumulative 
frequencies. Even in this crude model we can see how 
the distribution of the products is coming reasonably 
close to the distribution D(x). 


one digit | two digit infinite digit D(x) 
Galealie || fore. cum. can a: % % cum. % 
1871 asia eyo a/b 30.1 


SST 41.7 18.3 42.4 7.7 
4565 56.4 14.5 56.9 60.2 
5528 6883) le 7~< “68x6 69.9 
6309 18:0. 9.50 7eek 1.8 
6931 85a7. gal. Ong (eae 84.5 
7431 OV Gime lo Glas 90.3 
7811 96.6 4.6 96.4 95.4 
8100 100.1 3.4 99.8 100.0 
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Using a computer, we do the corresponding 8100 two 
decimal digit products to get the next two columns. 
Finally, using calculus, we get the infinite digit 
distribution in the next two columns. Comparing these 
results, we see that the two digit experiment gives 
results reasonably close to the infinite digit case so : 
that two digit arithmetic is not too severe an approximation! 
to actual practice. If we really doubted this, we could 
do the three digit case (using symmetry to save almost 
half the cases), but it does not seem worth it at this 
point. Finally, the last column gives the distribution 
D(x), and we see that starting with uniform distributions 
in x and y, the product x*y is already close to D(x). 


Suppose that instead of the uniform distribution g 
we took the factor x from the distribution D(x) and y from 
any other distribution you please. What would we get? 

To try it experimentally, we can take the midpoint number 
of each of the 90 intervals, 1.0 to 1.1, 1.1 to 1.2, 
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1.2 to 1.3,...,9.9 to 10.0 as a typical value for the 
: numbers in their corresponding intervals. The 
; probabilities of falling in these intervals are, of 


course, 
: Bite ( 
i ) 


D(i0.0) - D(9.9), 


for the x numbers, and the corresponding numbers for the 
y distribution whatever it is chosen to be. What will 
we get? For each y that you pick, the effect of 
multiplying those from the x distribution, because of 
the earlier result, 


D(kx) = D(k) + D(x) 


will merely reproduce the distribution D(x). Since 

this happens for each y value that you pick, regardless 
of the distribution of the y's, the result will be that 
the products will fall in the distribution D(x). Thus, 
if one factor comes from D(x) then the products will fall 
in the distribution D(x) regardless of what distribution 


the other factor comes from. This is called the ; ® 
"persistance of the distribution D(x)." If you try them, : 
computer experiments will bear this out. A particularly : 


spectacular illustration of this persistance is the 
example of forming a product of N factors 


X1*Xp*xX3* -e+ Xy 
If any one of the factors is from the distribution D(x), 
then the resulting product is also from the distribution. 


We now turn to division. First we show that if 
x is from Bin) then so is 1/x. We set x*y b (the 
number base). Then each of the following steps is 
obvious, after a little thought. Let P(x) be the 
distribution of 1/x. We have the definition 


D(x) = Prob{t < x} 
Then 
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Using this result we see that the behavior of 
mantissas in the process of division is reduced to the 
already solved case of their behavior in multiplication. 
Again, computer experiments, like those we did for 
multiplication, will verify this remark. When the 
original numbers x, y, and z, are all from the flat, 
equilikely distribution, the distribution of 


x*y/z 


is particularly close to the distribution D(x). 


Let us summarize what we have learned. First, 
the distribution D(x) that first arose when considering 
the distribution of physical constants also arises from 
the processes of multiplication and division provided 
we assume some reasonable distributions for the initial 
numbers. Thus, we have an alternate basis for expecting i 
to find the distribution D(x). Second, simple computer : 
experiments enable us to see how numbers combine; indeed, 
in some ways these experiments teach us more than do i 
the elaborate mathematical derivations, and in fact were 
the way the phenomenon was first discovered. 


EXERCISES 


1. Experimentally check that if x is taken from the 
distribution D(x), (meaning that the fraction of 
numbers less than x is (log x)/(log b), where by 
is the number base), then 1/x also has the same 
distribution. 


2. Extend the results of table 1 to three decimal digit 
arithmetic, and summarize as in table 2. 


3. Using the method of dividing the range into 90 
subintervals and selecting the mid-value as being 
typical of those in the subinterval, examine the 
distribution of the 8100 products x*y and the 
729,000 terms x*y/z assuming (a) each factor is 
selected from an equilikely distribution, (b) z is 
from D(z) and the others are from any other 
distributions you please. 


BY ROBERT TEAGUE 


Beginning a monthly column of discussion of 
programming languages. Professor Teague is 
at California State University, Northridge. 


The objectives of this column are: 


1) To have fun with the various programming 
languages. 


2) To encourage more communication among those 
persons interested in programming languages. 


3) To explore the facets and pitfalls of any and 
all programming languages. 


4) To promote the use of ANSI standards by 
uncovering the quirks due to the variations 
among compilers. 


Since this is the initial appearance of this column, it 

might be fun as well as educational to test our abilities 

in Algol, COBOL, Fortran, and PL/1 through a short quiz. 

The first person (determined by the earliest postmark) 

to submit correct answers to all questions will be : 
declared the winner and his name will be published in ; & 
a subsequent issue. Send your answers to: : 


Speaking of Languages... 
POPULAR COMPUTING 

Box 272 

Calabasas, CA 91302 


1. In PL/1, what does the following notation mean 
and 
what data type does it imply for K and A? : 


K—A 


2. In COBOL, what level number must be assigned to 
OVER-18? 


IF 6VER-18 THEN G@ T% DRAFT-STATUS. 


3. In Algol, what kind of a division is implied by = 
the following? E 


Fencervenasonsencons 


4, In Fortran, is the following statement legal? What 
type is implied for K? 


IF (K) G@ T@ 10 


5. Recode the following Algol statement into one PL/1 
assignment statement. 


L i= if A>B then 4 else 0; 


6. Recode the following Algol statement into one 
statement in Fortran, COBOL, and PL/1. 


switch A := L10, L20, 130, TAO, L50; 
goto A [I]; 


7. Recode the following Algol for statements into one : 
different statement in Algol, | COBOL, Fortran, and nd PL/1. 


a) for J :=M step 3 until 100 do (etc.) 
b) for L := 8, I++ while L<150 do (etc.) 


8. Recode the following PL/1 assignment statement into 
one statement in Algol, COBOL, and Fortran. 


A=B=4; 


9. Given the following PL/1 code, what will be the 
values of I and J after PL_1 is executed? 


DCL (A,B,C,D,E,I,J) FIXED BINARY; 


Nscspisove Coes we Se i = lp 
I = -A+B/C*(-D**2); 
PL1i: J=7(E&B<D); 


10. Given the following COBOL entries, what will be the 
values possessed by E(1,2,1), E(2,1,2), D(i,1), 
D(2,2), C(2)? 


02 FILLER PIC 99 VALUE 1. 
02 FILLER PIC 99 VALUE 2. 
O02 FILLER PIC 99 VALUE 3. 
: 02 FILLER PIC 99 VALUE 4. 
i 02 FILLER PIC 99 VALUE 5. 
i 02 FILIER PIC 99 VALUE 6. 
02 FILLER PIC 99 VALUE 7. 
02 FILLER PIC 99 VALUE 8. 

O01 B REDEFINES A, 

02 C OCCURS 2. 

03 D OCCURS 2. 
O4 E PIC 99 OCCURS 2. 
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The Raindrop Problem 


PROBLEM 


Given a square of side l. Three uniformly 
distributed random numbers are drawn. The first two 
are uniform in the range zero to one, and are taken as 
the X and Y coordinates of a point in the unit square. 
The third number is uniform in the range zero to 1/2, 
and is taken as the radius of a circle. How many such 
triplets must be chosen and plotted so that the unit 
square is completely covered? 


Flowcharting the procedure for selecting the 
required random numbers and describing (or plotting) 
the circles is relatively easy. The critical step is 
that of determining when the area covered by the circles 
also covers the complete square. 


There are many excellent algorithms for generating 
random numbers with a computer. Indeed, the literature 
of the field has only one topic that occurs more 
freyuently; namely, sorting techniques. 


If one has only a pocket calculator at his command, : 
and it functions only in floating decimal, many of the : 
standard algorithms are not applicable. For casual use : @ 
of random numbers, perhaps the old inner product method i 
is again of interest. As a computing algorithm, this 
method has been thoroughly discredited, for reasons 
given below, but as a hand method (where the user can 


monitor the output by eye) it has some merit. ® 


The method proceeds like this. Begin with two 
4-digit numbers. Form their product, which will have 
8 digits or less. Delete the low order two digits and 
retain the next four digits. Using this new number, 
and the one preceeding it, repeat the multiplication, 
obtaining the next 4-digit number in the sequence. For 
example: 


1234 
5678 
0066 
3747 
2473 
2663 
5855 
5918 
6498 and so on. 


As an algorithm to be used as a subroutine in a 
computer program, this scheme has two defects: 
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(1) The algorithm can degenerate to zeros at any 
time. The useful cycle length is thus unpredictable 
and is, of course, a function of the starting values. 
The cycle length to the point of zero degeneration is 
also usually relatively short. The degeneration point 
can be delayed by various tricks, but the method is 
inherently unstable. At least, however, the point at 
which the output goes to zero can be detected readily. 


(2) The scheme can degenerate in another way. 
Suppose that it has produced 10,000 good 4-digit numbers, 
but then the 10001ist and 10002nd numbers are the same as 
the 900lst and 9002nd numbers. From that point on, the 
same set of 1000 numbers will be generated over and over. 
This form of degeneration is difficult to detect, and 
vitiates the scheme as a suitable computer random number 
generator. 


Despite these serious defects, the scheme (which 
is the oldest aigorithm suggested as a random number 
generator) lends itself to hand calculation. The sample 
sequence started above, carried through 175 terms, 
develops the ten digits in this distribution: 


O ONAN FUWNPrO 
Ov 
=] 


i Testing these observed frequencies against the theoretical : 
i (70 in each case) gives a chi-squared value of 11.8. With: 
? 9 degrees of freedom, the probability level is around .25, : 
i which is acceptable for the frequency test. 


COMPUTIST: a specialist in computing. One who is 
: skilled in or practices the art of computing. One versed in 
computation. 
: COMPUTERIST: a specialist in computers. One skilled 
> and knowledgeable about the construction and design of 


: computers. 


Widespread use of desk calculators has revived interest 
in a lost art. Old tricks are suddenly new tricks again. : 
For example, on machines with sufficient accumulator capacity, 
sums of squares and cross products (of small numbers) can : 


be accumulated in one operation. Suppose we have this data: 
X ¥ 
a: >» 
4 5 
6 7 , 
8 9 i 


If the first pair is entered as a single number, 20000003, 
and that number is squared and summed, and the same procedure 
is followed for each pair, then the final result shows 


12000002800000164 


which is both sums of squares and twice the sum of the cross 
products. If this were done on an old rotary calculator, 
the secondary dial also shows 


200000024 


which gives the two simple sums. At present, this trick is 
not feasible on the electronic machines. 


On the other hand, elaborate tables that were built 
to facilitate desk calculator work are again valuable. One 
such table is 


Tables in Brief 

Mathematical Tables for Use with Desk 
Calculators 

Wayne White, June, 1951 

RM-660 

The RAND Corporation 

1700 Main Street, Santa Monica, CA 90406 

$2.00, 34 pages 


The booklet contains tables (most of them on one page <a 
of common logarithms, natural locarithms, exponentials, 10* 
sine, cosine, tangent (and their inverse functions) in both 
degrees and radians, hyperbolic functions, square roots, 
samma function, elliptic interrals, and Bessel functions-- 
and all of these arranged for ease of interpolation, both 
linear and quadratic. For example, the following are 
entries in the logarithm table: 


ie Log A B 
2.25 .35218252  .00193020 ~.00000428 
2.26 .35410844 .00192166 -.00000424 


To acquire log 2.2573, 


log (2.25+.0073) = log(2.25) + .73(.00193020 - 


3535884416 
( .3535893 true) 


Wt 


one calculates (for linear interpolation) ; 


.00000428 ) 


For quadratic interpolation: 


log(2.2573) = 


log(2.25) + AP + BP* 


where A and B are given in the table and P is the fraction .73. 


= .3535892852 
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